#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @File  : 坐标转换.py
# @Author: hukelin
# @Date  : 2019/7/18
# @Desc  :hello world!
 
'''
百度坐标转换
'''
import json
import urllib
import math
from urllib.request import urlopen, quote
import string
from math import pi,sin,cos
x_pi = 3.14159265358979324 * 3000.0 / 180.0
# pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 偏心率平方
# coordinate = []
# lng = []
# lat = []
# converted_lng = []
# converted_lat = []


def get_distance(point1,point2):
    lat1=point1[1]
    lat2=point2[1]
    lng1=point1[0]
    lng2=point2[0]

    EARTH_REDIUS = 6378.137

    def rad(d):
        return d * pi / 180.0

    radLat1 = rad(lat1)
    radLat2 = rad(lat2)
    lat_a = radLat1 - radLat2
    lng_b = rad(lng1) - rad(lng2)
    result_s = 2 * math.asin(math.sqrt(math.pow(sin(lat_a/2), 2) + cos(radLat1) * cos(radLat2) * math.pow(sin(lng_b/2), 2)))
    result_s = result_s * EARTH_REDIUS

    return result_s

 
def gcj02_to_bd09(lng, lat):
    """
    火星坐标系(GCJ-02)转百度坐标系(BD-09)
    谷歌、高德——>百度
    :param lng:火星坐标经度
    :param lat:火星坐标纬度
    :return:
    """
    z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi)
    theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi)
    bd_lng = z * math.cos(theta) + 0.0065
    bd_lat = z * math.sin(theta) + 0.006
    # print('gc==========>',lng,lat)
    # print('bd==========>',bd_lng,bd_lat)
    return {'lng':bd_lng, 'lat':bd_lat}
 
 
def bd09_to_gcj02(bd_lon, bd_lat):
    """
    百度坐标系(BD-09)转火星坐标系(GCJ-02)
    百度——>谷歌、高德
    :param bd_lat:百度坐标纬度
    :param bd_lon:百度坐标经度
    :return:转换后的坐标列表形式
    """
    x = bd_lon - 0.0065
    y = bd_lat - 0.006
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)
    gg_lng = z * math.cos(theta)
    gg_lat = z * math.sin(theta)
    return [gg_lng, gg_lat]
 
 
def wgs84_to_gcj02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
    """
    if out_of_china(lng, lat):  # 判断是否在国内
        return {'lng':lng, 'lat':lat}
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return {'lng':mglng, 'lat':mglat}
 
 
def gcj02_to_wgs84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    if out_of_china(lng, lat):
        return {'lng':lng, 'lat':lat}
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng

    return {'lng':lng * 2 - mglng, 'lat':lat * 2 - mglat}
 
 
def bd09_to_wgs84(bd_lon, bd_lat):
    lon, lat = bd09_to_gcj02(bd_lon, bd_lat)
    return gcj02_to_wgs84(lon, lat)
 
 
def wgs84_to_bd09(lon, lat):
    lonlat = wgs84_to_gcj02(lon, lat)
    # print('wgs84_to_bd09======>',lonlat)
    lg=lonlat['lng']
    lt=lonlat['lat']
    res=gcj02_to_bd09(lg, lt)
    # print(res)
    return res
 
 
def _transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
          0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret
 
 
def _transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
          0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret
 
 
def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    # print('out_of_china======>',lng,lat)
    return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)
 
 
# pathIn = r"你的文件路径\address_ln&la_Test.txt"
# pathOut = r"你的文件路径\result_ouput.txt"
# result = []
# count = 0
# if __name__ == '__main__':
#     with open(pathIn, 'r') as file:
#         lines = file.readlines()
#         for line in lines:
#             tmp = line.split(",")
#             lng.append(float(tmp[0]))
#             lat.append(float(tmp[1]))
#     try:
#         for i in range(len(lng)):
#             converted = bd09_to_wgs84(lng[i], lat[i])
#             with open(pathOut, 'a') as f:
#                 print(str(converted))
#                 f.writelines(str(converted) + "\n")
#                 count += 1
#             print("第" + str(count) + "条数据写入成功...")
#     except Exception as err:
#         print(err)
        
# ————————————————
# 版权声明：本文为CSDN博主「胡克林」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
# 原文链接：https://blog.csdn.net/weixin_42032107/article/details/97497327